function result = optim_svm(data, filename, job, options)
    if strcmpi(job, 'l1')
        L1 = createL1Program(data.X, data.y);
        L1.x = linprog(L1);
        d = floor(L1.prog_d/2);
        L1.w = L1.x(1:d) - L1.x(d+1:end-1);
        L1.b = L1.x(end);
        save(filename, 'L1');
        result = L1;
    elseif strcmpi(job, 'softl1')
        % example
        % l = optim_svm(data, 'result/res_margin2', 'softl1', struct('C',1))
        if length(options.C) == 1
            L1 = createSoftL1Program(data.X, data.y, options);
            L1.x = linprog(L1);
            [N,d] = size(data.X);
            L1.w = L1.x(1:d) - L1.x(d+1:2*d);
            L1.b = L1.x(2*d+1);
            L1.epsilon = L1.x(2*d+2:end);
            L1.C = options.C;
            model = L1;
            save(filename, 'model');
            result = L1;
        else
            for i = 1:length(options.C)
                C = options.C(i);
                L1 = createSoftL1Program(data.X, data.y, struct('C',C));
                L1.x = linprog(L1);
                [N,d] = size(data.X);
                L1.w = L1.x(1:d) - L1.x(d+1:2*d);
                L1.b = L1.x(2*d+1);
                L1.epsilon = L1.x(2*d+2:end);
                L1.C = C;
                model(i) = L1;
                
                display(sprintf('i=%d C=%f', i, C));
            end
            result = model;
            save(filename, 'model');
        end
    else
        result = 'not supported job';
    end
end

function p = createL1Program(X, y)
    [N, d] = size(X);
    for i = 1:N, X(i, :) = y(i)*X(i, :); end
    
    prog_d = 2*d+1;
    p.f = ones(prog_d, 1); p.f(prog_d) = 0;
    
    p.Aineq = zeros(N, prog_d);
    p.bineq = zeros(N, 1);
    for i = 1:N, 
        p.Aineq(i,:) = [-X(i,:) X(i,:) -y(i)];
        p.bineq(i) = -1;
    end
    
    p.Aeq = [];
    p.beq = [];
    
    p.lb = zeros(prog_d, 1); p.lb(prog_d) = -Inf;
    p.ub = Inf(prog_d, 1);
    
    p.x0 = zeros(prog_d, 1);
    p.options = optimset('Display', 'off')
    p.solver = 'linprog';
    p.prog_d = prog_d;
end

function p = createSoftL1Program(X, y, options)
    C = options.C;

    [N, d] = size(X);
    for i = 1:N, X(i, :) = y(i)*X(i, :); end
    
    prog_d = 2*d+1+N;
    p.f = ones(prog_d, 1); 
    p.f(2*d+1) = 0;
    p.f(2*d+2:end) = C/N;
    
    p.Aineq = sparse([-X X -y -eye(N)]);
    p.bineq = -ones(N, 1);
    
    p.Aeq = sparse([]);
    p.beq = [];
    
    p.lb = zeros(prog_d, 1); p.lb(2*d+1) = -Inf;
    p.ub = Inf(prog_d, 1);
    
    p.x0 = zeros(prog_d, 1);
    p.options = optimset('Display', 'off');
    p.solver = 'linprog';
    p.prog_d = prog_d;
end